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Abstract 

Current  methods  of  dim  object  detection  for  space  surveillance  make  use  of  a 
Gaussian  log-likelihood-ratio-test-based  detector  that  only  considers  the  probability  that  a 
single  point  in  an  image  resulted  from  the  propagation  of  a  point  source  through  space. 
This  paper  proposes  the  application  of  a  correlating  Gaussian  log-likelihood-ratio-test- 
based  detector,  which  makes  use  of  all  the  pixels  in  the  window,  to  the  arena  of  dim 
object  detection.  The  requirements  for  comparing  these  detectors  using  the  same 
threshold  of  detection  are  established.  The  performance  of  the  detectors  is  compared 
using  images  taken  by  the  Space  Surveillance  Telescope  in  New  Mexico  of  a 
geostationary  satellite  as  it  goes  into  and  out  of  eclipse  during  the  vernal  equinox  March 
13-15,  2012.  The  correlator  is  determined  to  offer  superior  perfonnance  based  on  a 
comparison  of  the  relative  probability  of  both  algorithms  to  detect  objects  near  the 
threshold  of  detection,  among  other  metrics.  The  additional  computational  cost  of 
implementation  of  the  correlating  detector  is  analyzed  and  detennined. 
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A  STUDY  OF  DIM  OBJECT  DETECTION  FOR  THE  SPACE  SURVEILLANCE 

TELESCOPE 


I.  Introduction 


1.  Problem  Statement 

Dim  object  detection  is  the  practice  of  deciding  whether  a  very  dim  object,  such  as  an 
asteroid  or  faint  satellite,  is  present  in  a  certain  location.  It,  by  necessity,  deals  with 
trying  to  detennine  the  presence  or  absence  of  low-level  signals  in  the  presence  of  a  high 
system  noise.  In  space  surveillance,  dim  object  detection  is  typically  applied  to  the 
challenge  of  locating  objects  that  can  be  approximated  by  point  sources,  such  as  space 
debris,  satellites,  asteroids,  and  other  near  earth  objects  (NEOs). 

This  research  deals  primarily  with  qualitatively  and  quantitatively  comparing  the 
efficacy  of  two  methods  for  dim  object  detection,  tenned  in  this  paper  the  point  detector 
and  the  correlator,  both  of  which  rely  on  the  likelihood  ratio  test  (LRT).  Both  methods 
are  a  form  of  Bayesian  Analysis,  as  they  rely  on  the  evaluating  whether  the  data  collected 
is  the  likely  result  of  a  pair  of  hypotheses:  one  in  which  an  object  is  present,  and  the 
other  in  which  an  object  is  absent.  The  point  detector  is  based  on  evaluating  the 
likelihood  that  a  single  pixel  centered  in  a  window  contains  an  object,  and  it  is  currently 
in  use  in  the  Lincoln  Near  Earth  Asteroid  Research  (LINEAR)  and  Space  Surveillance 
Telescope  (SST)  programs.  The  correlator  works  by  evaluating  the  likelihood  that  all  the 
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pixels  in  a  window  of  data  result  from  the  propagation  of  a  point  source  at  the  center 
pixel. 

2.  Problem  Importance 

Dim  object  detection  has  multiple  applications.  It  is  used  in  national  defense  for 
detecting  satellites.  It  is  used  to  detecting  space  debris,  which  threatens  both  civilian  and 
military  assets  of  both  the  United  States  of  America  and  other  countries.  These  are  two 
important  tasks  whose  utility  is  immediate.  Also,  it  is  used  in  meeting  NASA's  mandate 
to  detect  asteroids  [1],  which  has  use  both  in  scientific  research  and,  somewhat  less 
likely,  in  the  possible  mitigation  of  an  object  from  space  that  would  impact  Earth. 

Obviously,  detecting  dim  objects  based  on  simple  inspection  is  unrealistic;  if  the 
objects  were  easy  to  distinguish  then  they  wouldn’t  be  dim.  For  this  reason,  detecting 
dim  objects  is  best  done  with  digital  signal  processing-based  methods.  Asserting  that  an 
object  is  certainly  present  is  difficult;  even  more  so,  asserting  that  an  object  is  not  present 
on  a  given  trajectory  into  space  is  a  philosophical  conundrum.  Thus  statistically-based 
methods  that  decide  the  relative  likelihood  of  condition  are  appropriate. 

Given  the  expense  of  building,  maintaining,  and  operating  large  telescopes,  it 
behooves  those  trying  to  detect  dim  objects  to  implement  the  best  algorithm  available, 
since  limited  telescope  resources  are  otherwise  squandered.  To  that  end,  this  research  is 
concerned  with  making  the  best  use  of  already  existing  imagery  and  equipment  for  the 
purposes  of  dim  object  detection.  In  essence,  employing  a  better  detector  is  about  doing 
more  with  less,  which  is  important  during  this  time  of  fiscal  crises.  As  will  be  seen,  an 
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improved  detector  does  more  than  just  what  additional  telescope  assets  would  do:  it 
actually  allows  the  detection  of  objects  that  would  likely  not  be  detectable  using  optics 
currently  in  use.  For  a  given  probability  of  false  alarm,  a  better  detector  allows  for  a 
dimmer  object  to  be  detected  or  for  a  higher  likelihood  of  detecting  an  object  of  given 
magnitude  near  the  threshold  of  detection. 

3.  The  Space  Surveillance  Telescope  (SST) 

This  thesis  uses  imagery  of  objects  in  geosynchronous  orbit  that  was  taken  by  the 
Space  Surveillance  Telescope  (SST).  SST  was  developed  to  close  gaps  in  current  space 
surveillance  methods  by  doing  "broad-area  search,  detection,  and  tracking  of  small 
objects  in  deep  space[2]."  It  is  able  to  track  space  debris,  which  can  be  a  threat  to 
military  and  commercial  space  assets,  as  well  as  the  assets  themselves.  It  is  also  capable 
of  detecting  other  dim  objects,  such  as  asteroids,  and  so  it  is  able  to  help  fulfill  NASA's 
mandate  to  discover,  catalog,  and  track  90%  of  all  near  earth  object  s  (NEOs)  that  are 
larger  than  140  meters  in  diameter  [JJ. 

It  is  the  first  large  three-mirror  Mersenne-Schmidt  telescope.  Its  3.5  meter  primary 
mirror,  combined  with  the  f/1.0  design  allows  for  an  extremely  wide  field  of  view  that 
covers  an  area  approximately  three  times  larger  than  the  Ground-Based  Electro-Optical 
Deep  Space  Surveillance  system  (GEODSS),  which  is  the  best  existing  geosynchronous 
earth  orbit  sensor.  As  a  result  it  is  able  to  survey  a  large  amount  of  the  sky  at  once.  Its 
compact  design,  combined  with  a  highly  maneuverable  mount  and  precise  steering, 
allows  it  to  perfonn  a  broad  search  or  to  track  objects  as  they  travel  across  the  sky. 
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Figure  1.  View  of  SST’s  CCD  and  wide  field  of  view  camera  while  looking  down  the 
optical  path.  It  is  easy  to  distinguish  the  twelve  separate  sections  of  the  mosaic  focal 
surface  and  the  supporting  electronics  that  surround  it.  The  telescope  is  roughly 
horizontal  in  this  image,  so  the  rim  where  the  walls  of  the  house  meet  the  roof  is  visible 
to  both  the  left  and  right  of  the  telescope.  (http://www.darpa.mil/Our_Work/TTO/ 
Programs/Space_Surveillance_Telescope_(SST).aspx.) 


SST  employs  a  unique  focal  surface  charged-coupled  device  (CCD)  that  is  composed 
of  several  pieces  of  silicon  that  were  knit  together  into  one  mosaic  and  then  curved  to 
match  the  shape  of  the  wavefront.  Exposure  is  controlled  with  a  mechanical,  high-speed 
shutter,  rather  than  controlling  the  exposure  time  of  the  CCD  electronically,  in  order  to 
reduce  the  complexity  of  the  in-the-chip  electronics.  Although  well-executed,  the 
mechanical  limitations  of  the  high-speed  shutter  prohibit  taking  images  that  exhibit  a 
short-exposure  optical  transfer  function  (OTF)  of  the  atmosphere,  and  so  all  the  images 
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that  result  show  the  effects  of  an  average,  not  an  instantaneous,  atmosphere.  All  images 
used  in  this  research  were  taken  while  the  SST  was  located  at  Atom  Site  in  the  White 
Sands  Missile  Range,  above  8000  feet.  At  the  seeing  conditions  at  the  Atom  Site,  the 
CCD  is  sensitive  enough  to  be  able  to  detect  small  objects  in  geosynchronous  orbit, 
roughly  22,000  miles  away. 

SST’s  control  and  data  processing  system  is  able  to  process  the  images  in  near  real¬ 
time  when  using  the  Dirac  point  detector,  which  will  be  described  later.  The  telescope 
captures  up  to  1  terabyte  of  data  per  night[3].  As  will  be  seen,  the  research  that  is 
presented  here  can  be  run  on  the  current  data  processing  system.  Should  any  additional 
computing  resources  be  need  to  maintain  the  ability  to  run  in  real-time,  then  they  can  be 
added  in  parallel,  as  opposed  to  extending  the  depth  of  the  current  pipeline. 

4.  Past  Research 

This  section  provides  a  brief  review  of  literature  that  discusses  NEO  detection 
methods  and  differentiates  it  from  the  methods  presented  in  this  paper. 

a)  Detection  without  Foreknowledge  of  the  Background 

A  log-likelihood  ratio  test  used  by  the  Massachusetts  Institute  of  Technology  (MIT) 
Lincoln  Laboratory  in  the  Lincoln  Near  Earth  Asteroid  Research  (LINEAR)  program  is 
described  in  "Applying  Electro-Optical  Space  Surveillance  Technology  to  Asteroid 
Search  and  Detection:  The  Linear  Program  Results. "[4]  Some  SST  data  is  also 
processed  with  this  method.  That  likelihood  ratio  test  (LRT)  detector,  here-on  referred  to 
as  the  "point  detector,"  is  one  of  two  detectors  compared  in  the  research  effort 
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documented  herein.  Care  has  been  taken  to  faithfully  duplicate  the  detector  in  the 
manner  described  in  Pohlig's  writing,  which  is  discussed  in  the  next  section.  Particularly, 
the  same  threshold  of  detection  that  MIT’s  paper  describes,  which  determines  the 
probability  of  false  alarm,  is  employed,  since  this  allows  a  comparison  between  the 
proposed  correlator  and  currently  implemented  system.  While  LINEAR  depends  on  the 
presence  of  an  NEO  in  multiple  frames  in  order  to  find  asteroids,  the  method  being 
presented  in  this  thesis  is  strictly  single-frame  detection. 

b)  Detection  with  Foreknowledge  of  the  Background 

In  the  January  1989  issue  of  IEEE  Transaction  on  Aerospace  and  Electronic  Systems, 
"An  Algorithm  for  Detection  of  Moving  Optical  Targets"  is  presented  which  can  be 
applied  to  NEO  detection[5].  In  this  paper  Stephen  Pohlig  accurately  models  the 
interaction  of  photons  with  a  charge-coupled  device  (CCD)  as  a  collection  of  Poisson 
random  variables  [6].  However,  in  his  derivation  of  the  probability  of  false  alarm  he 
departs  from  a  strictly  Poisson  analysis  by  using  a  Gaussian  approximation  of  the 
normalized  log-likelihood  function.  The  method  relies  on  foreknowledge  of  the 
background,  which  does  not  make  it  applicable  to  use  with  SST  in  their  current 
processing  scheme. 

In  his  2012  master's  thesis  Near  Earth  Object  Detection  Using  a  Poisson  Statistical 
Model  for  Detection  on  Images  Modeled  from  the  Panoramic  Survey  Telescope  and 
Rapid  Response  System,  Capt  Curtis  Peterson  compares  the  performance  of  single-frame, 
log-likelihood  ratio  test,  single  point  detectors  that  are  based  on  either  Gaussian  or 
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Poisson  models  photon  arrival[7].  He  successfully  develops  a  true  Poisson  log-likelihood 
detector  in  which  he  shows  that  the  threshold  of  detection  can  be  assumed  and  then  tuned, 
rather  than  derived.  Capt  Peterson's  conclusions  suggest  that  further  exploration  of 
Poisson-based  models  will  be  very  fruitful.  The  methods  presented,  however,  were 
dependent  on  the  "master  sky  image"  that  the  Panoramic  Survey  Telescope  &  Rapid 
Response  System  (Pan-STARRS)  generates,  which  is  a  running  average  of  the  stellar 
background,  and  therefore  cannot  be  directly  applied  to  the  SST. 

5.  Preview 

The  remainder  of  this  document  will  be  structured  as  follows.  First,  a  model  of  the 
telescope  and  the  imagery  that  it  is  expected  to  produce  is  discussed.  Then  the 
experiment  process  and  rational  is  explained.  The  statistical  theory  that  underlays  the 
both  of  the  detectors  is  presented  before  going  into  the  details  of  how  each  of  the 
detectors  works.  The  logic  by  which  they  are  able  to  be  directly  compared  to  each  other 
is  then  presented.  The  results  of  applying  the  detectors  to  the  experimental  data  are 
presented,  with  an  emphasis  on  comparing  the  relative  probability  of  detecting  dim 
objects  near  the  threshold  of  detection  using  each  of  the  detectors.  The  computational 
cost  of  each  of  the  detectors  is  then  calculated  and  compared.  Conclusions  are  drawn 
from  the  results,  and  finally  limitations  of  the  research  are  discussed,  as  well  as  possible 
avenues  of  future  work. 
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II.  The  Model  and  the  Experiment 


1.  The  Model 

A  model  was  created  to  validate  the  images  taken  by  the  telescope  and  help  explain 
their  characteristics.  This  model  is  parameterized  to  approximate  characteristics  of  the 
Space  Surveillance  Telescope  (SST),  located  in  the  New  Mexico  desert.  It  is  an  FI 
telescope  with  a  diameter  of  3.5  meters  and  typical  Fried's  seeing  parameter  of  12  cm. 
The  atmosphere  is  considered  to  be  long  exposure  and  the  object  to  be  detected  is 
modeled  as  a  point  source. 

The  transfer  function  of  the  telescope  can  be  approximated  by  propagating  a  plane 
wave  through  a  binary  screen  that  represents  the  aperture  of  the  telescope  using  a 
Fraunhofer  approximation  and  the  Fourier  properties.  The  approximation  is  valid  since 
the  object  is  quite  distant  in  relation  to  the  aperture  of  the  telescope  such  that  the  distance 
from  any  point  to  the  aperture  to  the  object  can  be  considered  equidistant.  Thus,  the 
transfer  function  of  the  telescope  ,  Hopt(u,v) ,  is  given  by  Equation  (1)  [8] 


Hopt(u,v)  -  F 


\*{p(x>y)) 

|  2 

J7  r{p(x,y)}‘ 

> 

dudv 

(1) 


Here  (x,y)  are  indices  that  mark  off  p(x,y),  a  1320x1320  binary  screen  with  a  circular 

pupil  that  is  660  pixels  across  in  the  center  of  the  array.  F  is  the  Fourier  operator  for 

transforming  the  operand  into  ( u,v ),  the  frequency  space.  This  binary  screen  does  not 

account  for  any  obscurations  or  imperfections  in  the  telescope,  but  it  gives  a  fairly  good 

approximation  of  how  the  telescope  itself  affects  the  images,  as  can  be  seen  in  Figure  2. 
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Spatial  frequency  in  D/Xf 


Screen  Used  to  Model  Aperture 


Pixels 


Figure  2.  Transfer  function  of  the  telescope,  H  ,  (top)  that  results  from  the 

propagation  of  a  plane  wave  through  a  circular  aperture  (bottom)  using  the  Fraunhaufer 
approximation.  The  intensity  of  the  transfer  function  denotes  the  transmissivity  of  the 
telescope  to  different  frequencies  of  light,  in  terms  of  D  /  A/  . 
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The  transfer  function  of  the  atmosphere,  Hatm(vrad) ,  is  given  by  Equation  (2)  [6] 


HaJVrad)  =  eXP 


-3.44 

5/3' 

l  ro  2 

(2) 


where  1  is  the  wavelength  of  light  in  micrometers,  f  =  3.5  m  is  the  focal  length  of  the 
telescope  in  meters,  and  ro  =12  cm  is  Fried's  seeing  parameter  in  centimeters  (which  is 
typical  of  the  conditions  where  SST  is  located).  The  radial  frequency  vrad  is  defined  as: 

(3) 


Vrad  =  ^U2+V2 


The  resulting  pattern  is  graphed  shown  in  Figure  3. 


Figure  3.  Transfer  function  of  the  atmosphere,  Hatm ,  that  results  from  the  propagation  of 
a  plane  wave  through  the  long-exposure  atmosphere.  The  height  of  the  figure  denotes  the 
transmissivity  of  the  atmosphere  to  different  frequencies  of  light,  in  terms  of  D  /  Xf  . 
Note  that  the  scaling  is  different  from  the  other  transfer  functions  graphs. 
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Close  Up  of  OTF  of  a  Pixel,  H  , 

pixel 


Spatial  frequency  in  D/Xf 


Figure  4.  Transfer  function  of  the  atmosphere,  H.el  ,  that  results  from  the  propagation 

of  a  plane  wave  through  the  square  aperture  of  a  pixel.  The  intensity  of  the  figure 
denotes  the  transmissivity  of  the  pixel  to  different  frequencies  of  light,  in  tenns  of 
D/Af. 

Since  each  pixel  is  itself  a  square  aperture,  it  has  the  transfer  function  of  a  two- 
dimensional  rectangle  function  [8] : 

Hpixel  (u,v)  =  sinc(mv)  sinc(vvv)  (4) 

where  w  is  the  one-dimensional  width  of  the  square  pixel.  The  resulting  pattern  is  shown 
in  Figure  4. 

Since  the  system  is  linear,  the  transfer  functions  can  be  combined,  as  in  Equation  (5) 
H target ’v)  =  H atm(U’V)H opt(u,v)H pixel(u,v)  (5) 
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of  which  the  Fourier  transform  is 


Krge^y)  =  F  '  {HlargeM}  (6) 

However,  the  OTF  of  SST  is  undersampled  by  a  factor  of  30.  To  properly  model  this, 
h u,rge (X’  T)  's  skip-sampled  by  a  factor  of  30  to  an  1  lxl  1  array: 


h(x,y)  =  hlarge(x,y),  s.t.  x  =  30x,  y  =  30 y 
and  then  nonnalized,  so  that  the  PSF  of  the  models  is  shown  in  Equation  (8): 

h(x,y) 


Kodel^y)  ' 


y)dxdy 


(7) 


(8) 


OTF  of  the  Model,  H  , , 

model 


Spatial  liequency  in  D/A.I* 


Figure  5  Transfer  function  of  the  model,  Hmodel ,  that  results  from  the  propagation  of  a 

plane  wave  through  the  entire  system.  Note  the  heavy  pixelation,  which  results  from  the 
undersampling  of  the  image  plane.  The  intensity  of  the  figure  denotes  the  transmissivity 
of  the  model  to  different  frequencies  of  light,  in  terms  of  D  /  Xf  . 
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The  OTF  of  the  system,  then,  shown  in  Figure  5,  is  the  Fourier  transform  of  the  PSF: 

H model  (u,v)  =  F{ hmodel  (x,  V)}  (9) 

The  pristine  image  of  a  star,  ostar  (x,  y) ,  which  is  modeled  as  a  point  source,  is 

formed  by  placing  a  30  photon-count  pixel  in  the  center  of  a  11x11  zero  background. 
Propagating  this  image  through  a  telescope  or  atmosphere  is  mathematically  equivalent 
with  convolving  it  with  the  transfer  function  of  a  given  system,  and  the  result,  istar(x,y), 
is  found  by  [8] : 

i star  (x,y)  =  O  star  (x,y)* 

^ model  (x,y)  (10) 

which,  in  the  frequency  domain,  is  equivalent  to: 

1  star  («>V)  =  Ostar  (« .  V)H model  (W>  V)  (  1  1 ) 

Background  light  and  dark  current,  Ndark  (x,_y),  are  modeled  so  that  the  data, 
dstar  (x,  y) ,  has  a  Poisson  background  with  a  mean  of  ten.  The  measured  data  is  modeled 
by  using  this  sum  as  the  mean  of  a  Poisson  process  to  account  for  shot  noise: 

^  star  (x,y)  =  11  O  star  ( w ,  z)hmodel(w-  x,z-  y)dwdz  +  Ndark  ( x ,  y)  (12) 

One  thousand  instances  of  a  noisy  star  image,  dstar{x,y)  ,  were  created  in  this  way  to 

provide  a  sample  of  the  expected  images  that  SST  would  produce  of  dim  objects.  To 
provide  a  sample  of  a  background  noise  that  SST  would  see,  another  thousand  instances 
of  a  noisy  background  ddark  (x,  y)  were  created  in  the  same  manner  with  the  same  seed  to 

the  random  number  generator,  except  that  where  ostar(x,y)  contained  a  centered  point 
source  on  a  background  of  zeros,  odark  (x,  y)  was  entirely  zeros: 

ddark(x,y)  =  Ndark(x,y ) 
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(13) 


Example  of  ddark 


Pixels 


Example  of  d  t 

r  star 


Pixels 


1 star 


2  4  6  8  10 

Pixels 


Figure  6.  Examples  from  the  sample  model  data  sets  ddark  on  the  left  and  dstar  in  the 
center,  both  of  which  have  been  seeded  with  the  same  noise.  On  the  right  is  istar ,  the 
result  of  propagating  the  point  source  through  the  system  without  adding  any  noise.  The 
image  of  dstar  is  only  different  in  that  it  contains  istar .  The  intensity  of  the  point  source  is 

three  times  that  of  the  background.  All  three  images  have  been  scaled  equally  so  that 
their  relative  intensities  can  be  judged.  Significantly,  the  brightest  pixels  are  not  where 
the  star  is  located;  pixels  (1,8)  and  (11,11)  of  dstar  contain  the  brightest  pixels. 


Shown  in  Figure  6  is  an  example  from  both  sample  data  sets  which  has  the  same 
noisy  background;  that  is,  the  only  difference  is  that  the  example  of  dstar  has  dim  object 
in  the  center  pixel.  As  is  apparent  to  the  reader,  the  object  is  very  dim  in  relation  to  the 
background,  as  will  be  seen  as  typical  of  dim  objects  in  the  experimental  data  near  the 
threshold  of  detection. 

2.  The  Experiment 

To  perform  this  research,  data  was  needed  of  one  or  more  dim  objects  at  varying 
levels  of  intensity  so  that  the  detection  algorithms  could  be  compared.  It  was  necessary 
to  know  with  certainty  the  position  of  the  dim  object  that  the  detection  algorithm  was 
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being  used  to  detect  so  that  it  was  possible  to  differentiate  between  false  positives  and 
detections,  and  between  the  absence  of  an  object  and  a  failure  to  detect. 

Since  this  research  deals  with  the  cutting  edge  of  detection  ability  and  with  a 
telescope  that  is  capable  of  imaging  objects  that  have  never  been  seen  before,  using  a  star 
catalog  to  select  objects  was  deemed  insufficient,  both  because  of  the  difficulty  of 
locating  dim  objects  with  certainty  and  also  because  of  the  sheer  number  of  different 
objects  that  would  have  to  be  selected  and  located  in  order  to  characterize  the 
perfonnance  of  the  detectors  at  different  light  levels. 

Research  that  studied  the  glint  reflected  off  satellites  during  lunar  eclipses  [9]  inspired 
the  experiment,  and  a  glint  experiment  done  by  the  United  States  Naval  Observatory 
(USNO)  served  as  the  blueprint!  1  ()].  In  the  USNO's  experiment,  a  geostationary  satellite 
was  imaged  going  into  eclipse  during  the  autumnal  equinox.  Data  provided  by  the 
USNO,  which  is  shown  in  Table  1  and  graphed  in  Figure  8,  shows  how  a  rapid  series  of 
images  revealed  the  satellite  growing  dimmer  gradually  until  it  disappears  entirely. 

In  our  experiment  the  recently-constructed  SST  telescope  imaged  Anik  FI,  a 
geostationary  communications  satellite,  as  it  went  in  and  out  of  eclipse  for  eleven  nights 
between  28  February  and  23  March  2012  during  the  vernal  equinox.  The  satellite,  shown 
in  Figure  7,  orbits  at  107.3°  W  longitude  and  is  40.4  x  9  meters  at  its  largest  dimensions. 
The  satellite  was  actively  tracked  from  its  orbital  elements  throughout  the  entire  eclipse. 
The  result  is  a  dataset  in  which  an  object  of  decreasing  visual  magnitude  is  at  a  known 
location,  providing  the  opportunity  to  test  detection  algorithms  under  a  wide  variety  of 
light  conditions.  Since  the  satellite  was  tracked,  other  objects  in  the  sky  pass  through  the 


15 


Figure  7.  Artist's  depiction  of  Anik  F 1  in  orbit.  (Reproduced  without  permission  from 
http://www.boeing.com/defense-space/space/bss/factsheets/702/images_702fleet 
/0  lpr_0 1 507_hirez.jpg) 

frames  at  the  sidereal  rate.  This  situation  presents  the  observer  with  an  object  in  a 
known  location  that  has  an  intensity  that  varies  in  time,  which  allows  detection 
algorithms  to  be  tested  over  a  variety  of  light  levels.  These  attributes  make  this  type  of 
test  ideal  for  comparing  the  performance  of  different  detection  algorithms. 

Data  captured  once  per  second,  resulting  in  relatively  gradual  changes  between 
frames.  Both  binned  (2x2)  and  un-binned  pixel  data  was  collected,  but  only  binned  data 
is  included  in  this  analysis.  So  that  a  curved  imaging  surface  could  be  created  for  SST,  it 
is  compromised  of  several  separate  pieces  of  silicon  that  have  been  combined  together; 
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Table  1.  Digital  counts  for  Anik  FI  entering  the  terminator  during  autumnal  equinox. 
Provided  by  Dave  Monet,  USNO 


From  USNO  Telescope 
w/Orthogonal  Transfer  CCD* 

Frame 

UT 

Counts  above  sky 

2500 

25:31.1 
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2550 

40:41.1 
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6000 

Log(Counts)  =  228.6391  -  0.0091714*Time  (GMT) 
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10 

1000 
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Figure  8.  Data  points  and  approximate  light  curve  fit  of  the  digital  counts  for  Anik  FI 
entering  the  shadow  during  the  autumnal  equinox  12  October  2011. 
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typical  imaging  planes  are  flat  and  their  silicon  is  monolithic.  As  a  result,  each  of  the 
sections  of  the  imaging  plane  on  SST  have  slightly  different  responses,  and  there  are 
visible  seams  between  the  pieces.  In  order  to  accommodate  this,  data  from  only  one 
section  of  the  imaging  plane  was  used  in  processing  a  given  night. 

Data  from  the  nights  of  March  13-15,  2012  were  used  in  this  study,  which  correspond 
to  nights  073,  074,  and  075,  respectively.  These  three  nights  were  chosen  because  their 
data  is  the  most  comparable  to  each  other.  Binned  data  was  collected  on  these  nights. 
On  other  nights,  either  unbinned  data  was  collected  or  the  images  suffered  from  problems 
like  an  illuminated  moon.  On  the  first  two  nights  the  telescope  captured  the  satellite  not 
only  as  it  went  into  eclipse  but  as  it  came  out,  as  well. 

The  dataset  does  not  provide  a  basis  for  forming  receiver  operating  characteristic 
(ROC)  curves.  Not  only  is  the  amount  of  data  insufficient,  it  is  also  extremely  difficult  to 
assert  that  in  a  given  location  that  an  undetectably  dim  object  does  not  exist,  which  would 
be  a  necessary  condition  for  computing  the  probability  of  false  alarm.  Instead,  the 
probability  of  false  alarm  is  taken  to  be  fixed  and  determined  by  the  Gaussian  distribution 
of  the  log-likelihood  detector  [4]. 

III.  Theory  and  Methodology 

1.  The  Likelihood  Ratio  Test 

The  likelihood  ratio  test  (LRT)  detectors  used  in  this  research  are  binary  decision¬ 
making  processes  founded  in  Bayesian  analysis.  Use  of  the  LRT  for  detection  of  dim 
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objects  has  been  well  described  before[5][7],  and  the  details  of  this  section  deal  mainly 
with  their  implementation  in  this  particular  research.  LRTs  reach  a  conclusion  by 
comparing  the  relative  probabilities  that  the  resultant  image  was  the  product  of  two 
competing  hypotheses  to  a  threshold.  It  is  important  to  understand  that  the  detector 
cannot  respond  with  absolute  certainty  that  a  dim  object  is  or  is  not  present,  only  which 
scenario  is  more  likely.  The  threshold  is  used  to  tune  the  level  of  certainty  with  which 
the  detector  responds. 

The  satellite  or  dim  object  is  considered  to  be  present  in  a  pixel  d(x,y)  in  the  image 
when  the  likelihood  ratio,  A (x,y)  exceeds  a  threshold,  rj  ,  as  in: 


A(x,y) 


P[d(x,y)\Hl\  > 

P[d(x,y)\H0]< 77 

H0 


(14) 


The  two  conditional  probabilities  are  the  likelihoods  that  the  recorded  image  resulted 
from  the  hypotheses  that  the  satellite  is  present  (Ht)  or  not  (  Hg ). 

The  probabilities  are  detennined  based  on  the  statistical  characteristics  of  pixels  in  an 
NxN  window  of  the  image  that  is  centered  on  pixel  (x, y )  .  Pixels  in  this  window  are 
referred  to  by  the  ordered  pair  (w,  z) .  For  (x,  v)  to  be  perfectly  centered  in  the  window, 
it  is  necessary  that  N  e  odd  integers.  The  upper  and  lower  bounds  of  the  window  are 
defined  in  relation  to  (x,  y )  thus: 


whigh  =x  +  N/  2-1/2 

wioW  =  x~N/2  +  1/2 

zhiSh=y  +  N/2-1/2 

ziow=y~N/2+l/2 
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Since  N  e  odd,  the  addition  or  subtraction  of  1/2  is  necessary  to  ensure  that  the  bounds 
refer  to  integer-numbered  pixel.  It  can  be  verified  that  there  are  precisely  N  integers  in 
the  sets  that  span  [wlow,whi  h]  and  [zlow,zhi  h]  when  the  endpoints  are  included.  These 

variables  are  used  as  the  bounds  for  all  summations  and  products  over  a  window. 

As  was  discussed  in  the  previous  chapter,  the  noise  is  better  characterized  as  Poisson 
rather  than  Gaussian.  However,  Gaussian  distributions  well  approximate  Poisson 
distributions  with  means  greater  than  50,  and  this  research  focuses  on  a  comparison  using 
the  assumptions  that  have  already  been  made  in  the  currently  fielded  point  detector. 
Thus,  by  assuming  that  the  noise  in  the  image  is  Gaussian,  the  probabilities  that  the 
object  is  present  ( Ht )  or  not  ( H0 )  are  [7] : 


Vhigh  Zhigh 


P[d(x,y)\Hl]  =  n  n  -7^' " 

w=w,„„,  z=z,„„,  v  ZTCCJ 


d(w,z)-B(x,y)-0hj 


(16) 


and 


,vhigh  ^ high  1 

p[d(x,y)\H„\=  n 

w=w,  z=z.  V  Z*TC  (J 


2(7 ' 


(17) 


where  w  and  z  are  the  indices  that  number  pixels  in  a  square,  NxN  window  of  the  image 
centered  on  d (x,  y)  ,  B  is  the  expected  value  of  the  background  of  the  window,  a  is  the 
standard  deviation  of  the  background  of  the  window,  9  is  a  scale  factor  for  the  expected 
amount  of  light,  and  h  is  the  point  spread  function  of  the  system. 

Whereas  analysis  of  data  from  other  experiments  like  the  Panoramic  Survey 
Telescope  &  Rapid  Response  System  (Pan-STARRS),  which  is  able  to  make  use  of 
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foreknowledge  of  the  background,  the  background,  B(x,  y ) ,  is  approximated  in  this 


comparison  by  the  mean  of  all  the  pixels  in  the  window  centered  on  (x,  y) : 


B(x,y ) 


whigh  zhigh 

Z  Z 


^  ^low  Z  Zlow 


d(w,  z) 
N2 


(18) 


This  is  a  good  approximation  when  the  window  size  is  chosen  to  include  the  tails  of 
the  PSF  (point  spread  function)  so  that  the  intensity  of  most  of  the  pixels  in  the  window  is 
dominated  by  the  background.  Should  the  window  be  chosen  to  be  too  small,  then  the 
system  will  suffer  from  a  lack  of  contrast  between  the  dim  object  and  the  background. 
Should  the  window  be  chosen  to  be  too  large,  then  not  only  will  the  computational  cost 
increase  unnecessarily,  but  the  risk  increases  of  multiple  objects  being  included  in  the 
window.  The  inclusion  of  multiple  light  sources  in  the  window  results  in  a  higher 
computation  of  the  background  level,  resulting,  again,  in  a  decreased  contrast  and  a 
generally  skewed  decision-making  process.  The  window  used  to  process  all  the  data  in 
this  research  was  15x15  pixels. 

The  calculation  of  the  standard  deviation  of  the  window,  <r,  is  particular  to  the 
implementation  of  each  detector  and  will  be  defined  below  for  each  detector. 

Taking  the  natural  log  changes  the  products  to  summations,  removes  the  constants, 
and  results  in  a  lower  dynamic  range  for  the  function  rill.  Simplifying,  the  log- 
likelihood  ratio  test,  L  ( x,y ) ,  is: 


L  (x,  y)  =  In  A(x,  y)  = 
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z  z 

W=Wlaw  z=zhw 
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n 


(19) 
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2.  The  Threshold 


An  examination  of  the  units  in  Equation  (19)  shows  how  the  threshold  is  defined.  In 
the  numerator,  the  pixels  in  the  image,  d(w,z),  and  its  mean,  B ,  are  both  in  units  of 
photon  counts,  and  their  difference  represents  the  signal;  its  multiplication  by  the  transfer 
function,  h(w,  z ) ,  which  is  dimensionless  and  will  be  shown  to  sum  to  one,  does  not 
change  the  signal  fundamentally,  but  only  its  representation.  In  the  denominator,  the 
standard  deviation,  a(x,  y ) ,  is  also  in  units  of  photon  counts,  and  it  is  a  measure  of  the 
amount  of  noise  in  the  signal.  Thus  the  threshold,  77 ,  is  defined  in  tenns  of  signal  to 
noise  ratio  (SNR). 

Since  for  the  case  of  H0  either  detector,  lnA(x,y),  will  be  shown  to  be  a  zero  mean, 
unit  variance,  Gaussian  random  variable,  the  threshold  is  a  measure  of  how  likely  a  false 
positive  is,  in  terms  of  standard  deviations  away  from  the  mean  of  a  Gaussian 
distribution.  This  allows  the  performance  of  both  of  the  detectors  to  be  compared  on  the 
same  scale,  in  particular  when  77  =  6  .  For  this  threshold,  an  integration  of  the  Gaussian 
curve  shows  that  the  probability  of  false  alarm,  Pja ,  is 


e2dx  =  9.866e10 


(20) 


3.  The  Point  Detector 

The  "point  detector,"  Ls(x,y),  is  currently  in  use  by  LINEAR  and  SST,  and  the 
following  explanation  of  its  working  is  based  on  Stokes  description  of  it  [4].  It  is  the 
result  of  considering  only  whether  the  center  pixel  in  each  window,  (x,  y )  ,  resulted  from 
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propagation  of  the  dim  object.  This  is  equivalent  to  treating  the  transfer  function  of  the 


system  as  a  Dirac  delta  function,  5 ,  a  binary  screen  that  blocks  out  everything  except  the 
center  pixel,  (x,  y)  .  Accordingly,  the  sifting  property  of  the  Dirac  can  be  used  to  make 
the  following  simplification  for  the  point  detector: 


L  (r  V)=  2  ^  (d(w,z)~  B(x,y))hs(w+ x,z  +  y) 


w=wlow  Z=Zlow 


°s(x,y) 


(21) 


'g  J  Q d(w,z)-B(x,y))S(w  +  x,z+v )  =  (d(x,y)- B(x,y))  ^ 


W=Wlow  Z=Zlow 


°s(x,y) 
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For  the  point  detector,  <js  ,  the  standard  deviation  of  the  background,  is  given  by: 


°e(x,y)  = 
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(24) 


When  H0  holds,  the  structure  of  the  point  detector  is  easily  recognized  as  a  prime 
example  of  transforming  a  random  variable  into  standard  nonnal  distribution  via 
subtracting  by  the  mean  and  dividing  by  the  standard  deviation. 

To  verify  that  L^x,  y)  is  zero  mean,  consider  that  when  an  object  is  not  present  that 

the  expected  value  of  the  center  pixel  in  the  window  should  be  the  same  as  the  expected 
value  of  the  window  as  a  whole: 
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E[(d(x,y)\H0\  =  E[B(x,y)\H0]  =  B(x,  y ) 


(25) 


As  a  result,  the  expected  value  of  the  numerator  is 

E[{d{x,  y)  -  B{x,  y))\H0]  =  E  [B(x,  y)  -  B(x,  y)]  =  0  (26) 

To  show  that  the  expected  value  of  the  mean  is  zero,  it  is  sufficient  to  show  that  the 
expected  value  of  the  numerator  is  zero  [6].  Thus: 


E[Ls(x,y)\H0] 


E[(d(x,y)-B(x,y))\H0_ 

E[as(x,y) 


0 

E[as(x,y)] 


0 


(27) 


4.  The  Correlator 

In  contrast,  the  likelihood  that  all  pixels  in  the  window  resulted  from  the  propagation 
of  the  dim  object  can  be  included  by  using  the  PSF  as  the  transfer  function  of  the  system 
instead  of  the  Dirac.  This  detector,  Lco;r  (x,  y) ,  termed  the  "correlator,"  applies  the 
principles  of  correlation  to  general  technique  used  by  the  point  detector.  It  is  based  on  an 
LRT  that  considers  the  correlation  of  all  the  pixels  in  the  window  with  a  point  spread 
function  (PSF)  of  the  system,  hPSF ,  that  is  centered  on  the  window: 


L  corr(X,y) 
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W=Wbw  -=Zlow 


(d(w,z)-B(x,  y))hPSF  ( w  +  x,z+  v) 
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(28) 


In  this  research  the  PSF  was  approximated  in  a  simplified  manner.  For  each  night,  a 
faint  star  on  the  same  section  of  the  image  plane  mosaic  was  selected  as  a  point  source, 
and  its  image  is  used  as  the  basis  for  the  PSF.  Several  criteria  were  used  in  selecting  the 
star  so  that  the  resultant  PSF  was  a  good  approximation  of  the  transfer  function  of  the 


24 


system.  First,  the  star  had  to  be  of  an  intensity  such  that  the  response  of  the  charge- 
coupled  device  (CCD)  was  linear;  that  is,  it  had  to  be  bright  enough  to  overcome  most  of 
the  shot  noise  but  dim  enough  so  that  the  CCD  did  not  exhibit  any  saturation  or  charge 
overflow.  Second,  the  point  source  had  to  be  roughly  centered  on  a  single  pixel;  that  is, 
the  star  did  not  seem  to  fall  between  two  pixels  or  in  the  comer  of  a  pixel.  Third,  the  bulk 
of  the  power  had  to  be  confined  to  roughly  a  3x3  area;  in  this  way,  the  tails  of  the  PSF 
were  well  represented  in  the  15x15  window.  (Given  the  simulated  transfer  functions  of 
the  system  and  the  way  that  stars  appear  in  the  images,  an  object  that  appeared  to  be 
larger  could  likely  not  be  considered  a  point  source.) 

Noise  in  the  raw  image  of  the  star,  dPSF(x,y) ,  is  reduced  by  subtracting  min (dPSF(x,y)) , 

which  is  the  lowest  value  found  in  dPSP(x,y) ,  from  all  of  the  pixels  in  dPSF (x, y) .  This 
removes  as  much  of  the  background  as  possible  while  preserving  the  tails  of  the  PSF; 
since  the  shot  noise  that  dominates  the  image  is  strictly  additive,  doing  this  subtraction 
will  remove  only  noise  and  not  any  signal.  The  differenced  image, 
dPSF(x,  y)  -  min (dPSF(x,  y)) ,  is  then  normalized  by  dividing  by  the  double  summation  of 


25 


Figure  9.  The  columns  show  the  normalized  PSFs  used  for  each  of  the  three  nights  of 
data,  from  left  to  right:  073,  074,  and  075.  The  first  row  shows  the  pixelated  image.  The 
second  and  third  rows  show  the  cross  sections  of  the  PSFs  for  the  x  and  y  axis, 
respectively.  The  cross  sections,  which  are  taken  across  the  center  row  of  each  axis, 
show  that  the  bulk  of  the  energy  is  contained  in  a  3x3  square  in  the  center  of  the  window 
and  how  the  tails  roll  off  towards  the  edges  of  the  window. 


itself,  causing  the  summation  over  the  normalized  PSF  to  equal  one;  this  is  done  to 
maintain  unity  gain  in  the  system  so  that  multiplication  by  the  PSF  only  changes  the 
shape  and  does  not  amplify  or  attenuate  the  signal.  (Similarly,  it  can  be  seen  by 
inspection  that  the  sum  over  all  the  pixels  of  the  Dirac  function,  hs,  also  sums  to  one 
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since  only  the  center  pixel  of  fallows  transmittance  and  it  does  so  with  unity  gain.) 
Thus,  the  normalized  PSF,  hPSF(x,y) ,  is  computed  by 


hpsF(x,y)  —  ■ 


dPSF(x,y)-min(dPSF ) 

i,  N  _  _ 

Z  Z  ( d  PSF  (x,y)~  min(WreF  )  ) 


(29) 


Detailed  graphs  of  each  of  the  PSF's  are  shown  in  Figure  9. 

These  PSFs  are  taken  to  be  good  measurements  of  the  true  PSF  of  the  system,  rather 
than  estimates.  This  is  due  to  their  high  SNR,  which  results  from  the  stars  being  very 
bright  compared  to  the  background  noise. 

A  less  labor-intensive,  more  sophisticated  approach  of  detennining  the  PSF  could 
quite  possibly  yield  superior  results.  There  is  reason  to  suspect  that  since  the  telescope 
does  not  sample  the  data  at  the  Nyquist  rate  that  differences  in  the  aliasing  of  the  star 
used  as  the  PSF  and  of  the  dim  object  affect  the  perfonnance  of  the  correlator.  However, 
as  will  be  seen,  the  correlator  still  yields  appreciable  improvement  over  the  point  detector 
when  using  a  PSF  derived  from  the  method  described. 

As  previously  discussed,  a  comparison  of  the  two  detectors  while  using  the  same 
threshold  necessitates  that  both  detectors  need  to  be  zero-mean,  unit-variance  Gaussian 


random  variables  when  no  dim  object  is  present,  as  in  the  H0  hypothesis.  To  establish 
these  conditions  for  the  correlator,  recall  that  the  point  detector  is  defined  as 

|_  (Y  v)=  2  g  (d(w,z)-B(x,y))S(w  +  x,z  +  y)  _  (d(x,y)~ B{x,y))  ^ 


w~w!ow  z  zlow 


<*e(x>y) 
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and  that 


E\_ls(x,y)\H,t}  =  0 


(27) 


The  standard  deviation  of  the  correlator,  <jcorr  (x,  y) ,  will  be  defined  at  this  point  in 
terms  of  the  multiplication  of  the  standard  deviation  of  the  point  detector,  <Js{x,y),  and 
an  unknown  tenn,  crtest(x,y ) ,  which  will  be  solved  for  later: 

0,  y)  =  <*s  (x,  y)vtest  (X,  y )  (30) 

The  terms  for  the  background  and  the  standard  deviation  do  not  change  based  on  the 
indexing  of  the  variables  for  the  summation.  Holding  these  constant,  the  definition  of  the 
point  detector  can  then  be  encapsulated  with  the  definition  of  the  correlator: 


L  (x  y)=  2  'y  ( d(w,z)-B(x,y))hPSF(w  +  x,z  +  y ) 
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The  mean  of  the  correlator  can  now  be  determined  by  substituting  in  the  mean  of  the 
point  detector: 


E\Lcorr(x,y)\H0]  =  E 


^  L s  (w,  z)hPSF  ( w+x,z+y ) 
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Turning  attention  to  finding  the  variance,  the  second  moment  of  the  correlator  is 
calculated: 
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When  (Wj,Zi)  =  (w2,z2) ,  substituting  one  for  the  other  reduces  this  to 
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In  his  article,  Stokes  established  that  the  point  detector  is  unit-variance,  so  we  can 


substitute  in  £  L/(wl,z1)  |  i/0  =1  [4],  Thus: 
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Since  we  desire  for  Lcorr  (x,  y)  to  have  a  unit  variance,  we  define  it  as  such  and  solve 

for  <r*«(x>y)  : 
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By  using  substitution  the  variance  of  the  correlator,  acorr  (x,  y) ,  can  be  fully  defined 
as: 

°corr  (x,  y )  =  as  (x,  y)crtest  (x,  y )  (30) 

whiSh  -ugk 

Zh  X  h2PSF(w+x,z  +  y)  (42) 

W=Wjow  Z=Z/ow 

Thus  the  both  the  point  detector  and  the  correlator  are  established  as  a  zero-mean, 
unit-variance  Gaussian  random  variable  and  can  be  compared  to  the  same  threshold  on 
equal  footing. 

5.  Search  Methodology 

The  image  of  the  satellite  moves  perceptibly  from  frame  to  frame.  Therefore,  the 
detectors  are  tried  on  all  the  pixels  in  a  6x6  search  grid  that  is  roughly  centered  on  the 
satellite.  The  satellite  is  considered  to  be  detected  when  it  is  found  in  any  pixel  in  the 
search  grid,  which  is  the  equivalent  of  applying  a  logical  “or.”  When  the  satellite  is 
detected  in  multiple  pixels  in  a  single  frame,  the  SNR  and  probability  of  detection 
calculations  shown  later  are  based  on  the  pixel  in  the  frame  that  had  the  highest 
probability  of  detection. 

For  the  purposes  of  this  analysis,  the  eclipse  is  considered  to  last  from  the  first  non- 
saturated  frame  to  the  last  non-saturated  frame.  The  algorithm  was  only  used  to  detect 
objects  during  this  defined  eclipse  period  for  a  two  reasons.  First,  the  objects  are  easily 
detectable  in  saturated  frames  by  simply  searching  for  pixels  that  have  an  intensity  that 
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equals  the  maximum  value  of  the  ADC,  so  it  is  a  trivial  case.  Second,  at  saturation 
charge  spillover  into  nearby  pixels  distorts  the  data  and  saturation  is  in  itself  a  significant 
non-linearity,  both  of  which  challenge  the  assumption  that  the  system  as  a  whole  is 
essentially  linear. 


IV.  Results 


1.  Detections 

The  correlator  showed  an  appreciable  gain  in  the  number  of  detections  over  the 
region  in  which  the  two  detectors  cross  the  threshold.  Even  though  the  experiment 
successfully  showed  the  satellite  pass  through  a  large  range  of  illumination,  the  time  in 
which  it  passed  through  the  threshold  was  relatively  small  compared  to  the  size  of  the 
data  set,  and  so  there  are  a  limited  number  of  frames  in  which  a  two  detectors  gave  a 
different  answer.  The  output  of  each  detector  as  the  satellite  goes  into  eclipse  is  plotted 
in  Figure  10.  The  correlator  did  not  universally  perform  better,  but  an  examination  of  the 
detections  near  the  threshold  shows  that  the  correlator  generally  detected  the  satellite 
deeper  into  the  eclipse.  Though  instructive,  simply  noticing  the  additional  frames  of 
detection  does  not  convey  much  knowledge  about  the  relative  perfonnance  of  the 
detectors.  The  small  number  of  frames  during  which  the  detectors  straddle  opposite  sides 
of  the  threshold  makes  it  difficult  to  draw  well-supported,  quantifiable  conclusions. 
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Figure  10.  Detections  as  the  satellite  goes  into  eclipse  for  all  three  nights  for  both 
detectors  near  the  crossing  of  SNR  6.  The  left  and  right  columns  show  the  results  of  the 
point  detector  and  correlator,  respectively.  From  top  to  bottom  the  three  nights  of  data 
are  073,  074,  and  075.  The  number  of  additional  frames  that  the  correlator  was  able  to 
detect  the  satellite  in  may  seem  insignificant,  but  the  span  of  frames  shown  represents  a 
significant  drop  in  apparent  magnitude. 


32 


2.  Comparison  of  Signal  to  Noise  Ratio 


As  has  been  established  in  the  discussion  of  the  threshold,  the  output  of  the  detectors 
for  a  given  pixel  directly  corresponds  with  the  signal  to  noise  ratio  (SNR)  for  that  pixel. 
The  algorithms  register  a  detection  when  the  SNR  of  a  pixel  rises  above  the  threshold. 
Figure  1 1  shows  the  SNR  of  both  detectors  as  the  satellite  goes  into  eclipse  on  night  073. 
The  point  detector  can  occasionally  be  seen  to  exhibit  a  higher  SNR  than  the  correlator, 


a 


Figure  11.  Signal  to  noise  ratio  for  both  detectors  as  the  satellite  goes  into  eclipse  on 
night  073.  The  correlator  is  depicted  in  red  and  the  point  detector  is  depicted  in  blue. 
The  horizontal  green  line  shows  the  threshold  of  SNR  6.  The  top  graph  (a)  shows  the  full 
data  set  and  the  bottom  graph  (b)  shows  a  close  up  as  the  satellite  fades  beneath  the 
threshold  of  detection.  The  sharp  spikes  in  (b)  around  frames  900  and  950  are  due  to  star 
crossings. 
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but,  on  average,  the  correlator  comes  out  on  top.  The  correlator  also  seems  to  experience 
less  perturbation  when  other  sources  of  light  interfere. 

The  discussion  that  now  follows  concentrates  on  the  region  of  the  datasets  where  the 
output  of  the  detectors  is  near  the  threshold,  since  those  data  points  are  the  ones  that 
allow  for  differentiation  of  the  detector's  perfonnance.  The  reason  is  that  objects  that  are 
much  brighter  are  not  interesting  since  they  are  easily  detectable  by  both  algorithms,  and 
objects  that  are  much  dimmer  are  not  interesting  since  the  probability  of  detection  by 
either  algorithm  approaches  zero. 

3.  Relative  Probability  of  Detection 

Attempting  to  compute  the  probability  of  the  detection  directly  poses  several  major 
hurdles.  Most  significantly,  it  requires  foreknowledge  of  the  presence  of  an  object  that  is 
difficult  to  detect.  The  light-gathering  power  of  the  Space  Surveillance  Telescope  (SST) 
makes  it  possible  to  see  objects  that  have  not  previously  been  detected.  Were  this 
research  conducted  using  images  collected  by  a  small  aperture  telescope,  it  might  be 
possible  to  rely  on  a  catalog  to  indicate  where  stars  are  that  would  appear  dim  to  that 
telescope;  this  is  obviously  not  a  possibility  for  SST. 

Similarly,  attempting  to  compute  the  probability  of  detection  directly  is  also 
extremely  difficult.  To  do  so  requires  absolute  knowledge  that  there  is  absolutely  nothing 
in  a  given  window:  a  philosophical  conundrum.  Instead  a  threshold  of  detection  is 
chosen  in  a  way  that  detennines  the  probability  of  false  alann  based  on  assumptions  in 
the  Gaussian  model. 
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None  of  the  differences  between  the  point  detector  or  correlator  change  the 
probability  of  false  alarm.  Holding  this  variable  of  perfonnance  constant  allows  for  the 
algorithms  to  be  compared  with  relative  ease.  While  the  absolute  probability  of  detection 
cannot  easily  be  detennined  for  either  method,  it  is  possible  to  determine  the  probability 
of  detection  of  one  as  a  function  of  the  probability  of  detection  of  the  other  using  an 
analysis  of  the  signal  to  noise  ratios  (SNRs). 

The  starting  frame  for  the  graphs  that  will  be  shown  in  this  section  is  when 
probability  of  detection  for  the  two  detectors  starts  to  recede  from  one,  and  the  ending 
frame  is  when  the  probability  of  detection  for  the  two  detectors  approaches  zero. 
Including  points  outside  of  this  range  results  in  clutter  near  the  endpoints  that  obscures 
the  features  of  the  graphs  and  does  not  add  any  meaningful  infonnation. 

For  a  given  object  with  a  specific  SNR,  the  probability  density  function  for  the 
probability  of  detection  is  predicted  by  a  unit-variance  (cr),  Gaussian  distribution, 
/(x;//,cr2): 


/(x;//,cr)  =  - 


1 


If  X~jU 

2v  a 


(43) 


crV27r 

Since  the  Gaussian  distribution  is  symmetric,  when  the  measured  SNR,  x ,  is  set  to 
equal  the  mean,  ju  ,  the  probability  density  reaches  a  maximum  value. 

Obviously,  the  likelihood  of  detection  increases  as  the  threshold  SNR  gets  lower,  just 
as  the  probability  of  false  alarm  also  increases.  In  the  same  way,  the  probability  of 
detection,  Pd ,  decreases  as  the  threshold  SNR  is  raised.  Thus,  the  PDF  should  be 
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integrated  beginning  from  a  threshold  SNR  of  positive  infinity,  towards  a  lower  bound 
threshold,  r  ,  as  in: 

Pd- j  f(x',jU,a2)dx  (44) 

Checking  the  upper  and  lower  bounds  for  the  threshold  shows  that  this  is  sensible. 
When  the  threshold  of  detection  is  set  impossibly  high,  such  that  r  =  go  ,  we  see  that 

(•go  0 

Pd  =  [  /  (x;/u,  cr"  )dx  =  0  (45) 

J  T= 00 

If  t  —  —00  then  the  SNR  is  set  as  low  as  possible  so  that  there  is  no  requirement  for 
any  signal  at  all;  in  this  case  the  integration  is  defined  as  one: 

Pd=\  f{x',ju,(7~)dx  =  \  (46) 

Jr=— 00 

Taking  notice  that  the  above  integration  are  not  from  "left  to  right,"  the  notation  can 
now  be  simplified  by  using  the  definition  of  the  cumulative  distribution  function  (CDF) 
of  the  Gaussian  distribution,  F(r;  //,cr ) ,  and  subtracting  it  from  one  [121: 

Pd  =  1  -F(t;  /j,  cr2 )  (47) 

Recall  that  the  mean  of  a  detection  algorithm  is  the  threshold  at  which  an  object  with 
a  specific  SNR  is  equally  likely  to  be  detected  as  not.  The  relative  probability  of 
detection  of  the  two  algorithms  can  then  be  compared  for  the  same  threshold,  r  ,  since  the 
sufficient  statistics  of  the  two  algorithms  are  unit  variance  Gaussian  random  variables  and 
they  are  being  compared  in  detecting  the  same  object;  the  only  difference  is  the  mean  of 
the  distributions. 
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Consider  a  single  data  point  in  which  the  threshold  is  SNR  five  and  the  mean  of  the 
correlator  is  SNR  six  for  an  object  with  SNR  six.  The  relative  probability  of  detection  for 
the  point  detector  is: 

Pds  =l-F(r;//j,c2)  =  l-F(6,5,l)  *1-0.83  =  0.17  (48) 

Likewise,  the  relative  probability  of  detection  of  the  correlator  is: 

Pdcorr  =  l ~  F(r;  <ucorr ,  a2 )  =  1  -  F(6, 6, 1)  =  1  -  0.5  =  0.5  (49) 

Each  frame  presents  both  algorithms  with  the  same  snapshot  of  an  object.  To 
simplify  the  computations,  the  threshold  of  comparison  is  chosen  as  six.  As  the  Gaussian 
distribution  swells  at  the  mean,  it  is  assumed  that  the  mean  of  each  algorithm  is  the  SNR 
at  which  the  object  was  detected. 

Having  covered  the  general  theory,  the  notation  is  now  expanded  to  include  an  index 
to  denote  the  frame  from  which  the  SNR  is  derived  and  for  which  the  probability  of 
detection  is  being  computed.  Then,  for  frame  f  the  probability  of  detection  for  each 
detector  is 

Pds  (f )  =  1  -  F(x(  f );  jus  (f ),  rr2 )  =  1  -  F(6,  SNRs  (f ),  1)  (50) 

and 

Pdcon.  (f )  =  1  -  F(x(  f );  Mcon.  (f ),  a2 )  =  1  -  F(6,  SNRorr  (f ),  1)  (51) 

This  comparison  only  yields  a  relative  comparison  of  the  probability  of  detection  of 

each  algorithm.  Since  it  depends  on  the  threshold  of  detection  being  the  same  for  both 

algorithms,  and  since  both  algorithms  typically  detect  objects  that  are  much  brighter  than 

the  threshold,  it  gives  an  estimate  of  the  relative  increase  of  the  objects  that  will  be 

detected  around  SNR  six,  not  an  increase  in  the  objects  that  will  be  detected  in  total.  This 
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is  reasonable,  since  we  would  not  expect  to  suddenly  be  able  to  detect  a  larger  number  of 
relatively  bright  objects  by  changing  to  a  different  algorithm. 

Since  integration  of  the  Gaussian  curve  is  very  sensitive  to  the  placement  of  the 
mean,  computing  the  probability  of  detection  in  this  manner  across  a  region  of  the  dataset 
yields  a  very  noisy  signal.  Implementing  a  low-pass  filter  on  the  SNR  signal  that  before 
it  is  fed  into  the  integration  removes  the  high-frequency  noise  and  gives  a  better  picture 

of  the  perfonnance  trends  of  the  two  detectors.  The  low-passed  signal,  SNR(f ) ,  is 
generated  by  taking  the  running  average  of  15  frames  of  dataf  131: 

I*(g) 

Sm(£)  =  &— -  (52) 

The  affect  of  this  low  pass  filter  on  the  SNR  is  shown  in  Figure  12,  which  shows  the 
instantaneous  (top)  and  low-passed  (bottom)  SNR's  for  both  detectors  as  the  satellite  goes 
into  eclipse  on  night  073. 

Low-passed  versions  of  the  probability  of  detection  are  then  defined  based  on  the 
low-passed  SNR(f ) : 


Pds(f)  =  l-F(6,SNRs(f),l) 

(53) 

Pdcorr(f)  =  l-F(6,SNRcorr(f),l) 

(54) 

Possibly  the  most  illuminating  comparison  is  when  the  probability  of  detection  for 

each  of  the  algorithms  is  compared  per  frame.  Figure  13  shows  plots  of  Pds  and  Pdcorr 

over  a  span  of  adjacent  frames  as  the  probability  of  detection  sweeps  from  near  zero, 

through  the  threshold  of  detection,  to  practical  certainty  for  the  two  algorithms  as  the 
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satellite  goes  into  and  out  of  eclipse.  Also  shown  is  a  comparison  of  the  low-passed  Pds 


and  Pdcorr  ,  which  give  a  clearer  picture  of  the  generalized  performance  of  the  correlator 

on  a  given  eclipse  event.  Since  the  point  detector  is  currently  in  use  at  LINEAR  and  SST 
it  is  shown  as  the  baseline  for  performance. 


a 


b 


Frame 


Figure  12.  Instantaneous  (a)  and  low-passed  (b)  SNR  of  the  detectors  as  the  satellite 
goes  into  eclipse  on  night  073.  The  correlator  is  depicted  in  red  and  the  point  detector  is 
depicted  in  blue.  The  horizontal  green  line  shows  the  threshold  of  SNR  6.  Note  how  the 
low-passed  SNR  reveals  an  undulating  pattern  of  lower  sensitivity  that  affects  both 
detectors. 


39 


Pd  of  the  Compared  Detector  Pd  of  the  Compared  Detector  Pd  of  the  Compared  Detector 


073  In  Instantaneous  Pd-Pd  Curve 


073  In  Low-Passed  Pd-Pd  Curve 


073  Out  Instantaneous  Pd-Pd  Curve 


073  Out  Low-Passed  Pd-Pd  Curve 


074  In  Instantaneous  Pd-Pd  Curve 


o  °  o 


o  o 
)  o  °  o  , 


°  o  o 
o  o  a 

o  0  o£° 

°°  y 

& 


.2  .4  .6  .8  1 

Pd  of  the  Point  Detector 


o  .8 


(X 

g 

c 

U 


074  In  Instantaneous  Pd-Pd  Curve 

„  °  o  „  o«©G 


o  0  o 
o 

8 


o°  8 


°° 


.2  .4  .6  .8  1 

Pd  of  the  Point  Detector 


40 


o  •* 
u 

Q 

13 

C3 

B. 

g 

o 

U  .4 

4> 

£ 

o 

2  .2 


Of 


074  Out  Instantaneous  Pd-Pd  Curve 

°  °  °  o°  °  °  '-■**9 

3_  „  o 


<P 


.2  .4  .6  .8  1 

Pd  of  the  Point  Detector 


075  In  Instantaneous  Pd-Pd  Curve 


o  .8 

£ 

l-6 

g 

^  .4 


■y 


✓ 


<• 

> 


"0  .2  .4  .6 

Pd  of  the  Point  Detector 


Q 

r3 


o 

'TD 

Oh 


074  Out  Low-Passed  Pd-Pd  Curve 

Cfc>  O  OO  O  O  O  <£"3D» 

°  7^ 

o  <P 


,  <p 


’  .4  .6  .8  1 

Pd  of  the  Point  Detector 


075  In  Low-Passed  Pd-Pd  Curve 


5  .8 

£ 

•8  .6 


g 

5 

0  .4 


Jr 


!  .4  .6 

Pd  of  the  Point  Detector 


Figure  13.  Comparison  of  the  probability  of  detection  for  both  detectors  on  all  three 
nights.  Graphs  are  labeled  "In"  or  "Out"  for  whether  the  satellite  is  going  into  or  out  of 
eclipse,  respectively.  In  the  left  column  the  instantaneous  probabilities  of  detection  are 
Pds  in  blue  and  Pdcorr  in  red.  In  the  right  column  the  low-passed  probabilities  of 

detection  are  Pds  in  blue  and  Pdcon.  in  red,  both  of  which  are  the  product  of  low-pass- 

filtered  data  that  uses  a  15  frame  running  average.  The  data  points  used  for  these  graphs 
are  those  which  are  near  the  threshold  of  detection  for  each  of  the  algorithms;  they  were 
arranged  in  order  of  increasing  probability  of  detection  when  computed  by  the  point 
detector. 
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Figure  14.  Co-plotted  curves  of  Pdcorr  for  all  five  recorded  eclipse  events  over  nights 
073,  074,  &  075. 


A  condensed  view  of  the  relative  probability  of  detection  of  both  detectors  over  all 
five  eclipse  events  over  the  three  nights  is  shown  in  Figure  14.  The  differences  in  the 
low-passed  curves  shows  that  the  performance  gains  of  the  correlator  vary  significantly 
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from  a  small,  though  still  meaningful,  improvement  as  the  satellite  enters  the  eclipse  on 
night  073  to  over  a  near-guaranteed  detection  for  the  correlator  while  the  point  detector  is 
still  struggling  to  detect  the  satellite  as  it  exits  eclipse  on  the  night  of  074.  Over  the  range 
in  which  the  point  detector  has  between  a  20%  and  a  50%  chance  of  detecting  the 
satellite,  the  low-passed  probability  of  detection  for  the  correlator  ranges  between  just  as 
likely  to  detect  to  four-and-a-half  times  as  likely  to  detect.  Generalizing  the  five  curves, 
the  correlator  averages  a  75%  chance  of  detection  when  the  point  detector  has  a  50% 
chance  of  detection,  and  the  correlator  averages  a  50%  chance  of  detection  when  the 
point  detection  has  a  20%  chance  of  detection. 

4.  Cost  to  Implement 

The  computational  burden  of  each  of  the  algorithms  is  computed  by  counting  the 
number  of  additions  and  multiplications  that  have  to  be  perfonned  by  the  algorithms  on 
average.  While  certain  programming  techniques  can  be  result  in  efficiencies  for 
perfonning  addition  and  subtraction,  modern  superscalar  processors  perform  floating 
point  multiplication  as  easily  as  floating  point  addition  and  subtraction.  Incredibly,  "the 
complexity  of  square  root  evaluation  or  division  is  the  same  as  that  of  multiplication" 
when  the  Newton-Raphson  method  is  used  1 141.  Thus  floating  point  operations  for 
addition,  subtraction,  multiplication,  division,  and  square  roots  are  assumed  to  impose  the 
same  computational  burden.  Furthermore,  it  is  assumed  that  element-wise 
multiplications  of  matrices  are  performed  using  a  matrix  or  vector  processing  unit  in 
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which  loading  the  appropriate  data  based  on  a  region  does  not  impose  an  appreciable 
computational  burden. 

First,  consider  the  formula  for  the  point  detector,  L s(x,y) : 


(d(x,y)-B(x,y)) 

<rs(x,y) 


(22) 


where  (x,y)  is  the  pixel  being  examined,  d(x,  y)  is  the  data  point  in  that  pixel,  B(x,y)  is 
the  background  in  an  NxN  window,  and  crs(x,y)  is  the  standard  deviation  of  the  data  in 
the  window.  Given  B(x,y)  and  crs(x,y),  evaluating  the  point  detector  requires  only  a 
subtraction  and  division,  which  is  two  operations. 

The  background,  B(x,y) ,  is  approximated  for  each  window  by  taking  the  mean: 

wHSh  zhigh 

X  X  d(w’z) 

B(x,y)  =  —-^2 -  (18) 


Evaluating  the  numerator  requires  N2  - 1  additions,  evaluating  the  denominator  requires 
a  single  multiplication,  and  the  taking  the  ratio  requires  a  single  division.  As  N 
increases,  the  number  of  operations  needed  in  determining  it  rapidly  converges  to 

(a2-1)  +  1  +  1«  N2  (55) 

Likewise,  a ,  the  standard  deviation  of  the  background,  is  given  by: 


whigh  zhigh 

Z  Z  d2(w,z ) 


z~zlow 

N2 


B(x,yf 


(24) 
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Squaring  d(w,z )  for  each  pixel  requires  A2  multiplications,  and  the  summation  itself 


requires  N2  -1  additions.  Dividing  by  A 2  adds  only  two  operations,  as  with  evaluating 
the  background.  Subtracting  the  background  requires  only  a  single  operation  since  it  has 
already  been  computed,  and  one  more  again  for  the  square  root.  As  A  increases,  the 
number  of  operations  needed  in  determining  a  rapidly  converges  to 

N2  +  (A2  —  l)  +  2  +  l  +  l«  2N2  (56) 

Recalling  the  definition  of  the  point  detector,  altogether  evaluating  the  point  detector 
requires  N2  for  the  background,  2N2  for  the  standard  deviation,  and  an  additional 
subtraction  and  multiplication,  so  the  total  number  of  operations  needed  converges  to 

N2  +  2N1  +  2  «  3  A2  (57) 

Second,  consider  the  formula  for  the  correlator: 


|_  (x  y)=  '2  g  ( d(w,z)-B(x,y))hPSF(w  +  x,z  +  y ) 


W=Wlow  Z=Zlow 


_ PSF ' 

^corr(^y) 


(28) 


where  ( x,y )  is  the  center  pixel  of  the  Ax  A  window  being  correlated,  w  and  z  are 
indexes  in  the  window,  d{w,z)  is  the  data  point  at  a  specific  index,  B(x,y)  is  the 
background  in  the  window,  hPSF(w+ x,z  +  v) refers  to  the  pixel  of  the  PSF  for  a  point 
source  of  the  system  that  is  centered  at  (x,y)  and  has  been  nonnalized  for  the  window, 
and  (7corr  (x,  y)  is  the  standard  deviation  of  the  window. 

From  Equation  (30),  the  standard  deviation  for  the  correlator  is  found  by  multiplying 
the  standard  deviation  for  the  point  detector  by  <Jtest  (x,  y)  ,  which  is  defined  as: 
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^(x,y) 


(41) 


Wlow  Zhigh 


X  Tj  h2psF(wi+x’zi+y) 


W\  =Wlow  Z1  =Zlo\ 


Squaring  h2 PSF  (wl  +x,z1  +  y)  element-wise  and  taking  the  summations,  followed  by 
the  square  root,  requires  2N  operations.  However,  since  Equation  (41)  is  not  dependent 
on  d(w,z )  or  any  tenns  that  are  dependent  on  d(w,z),  it  does  not  change  as  (x,y) 
changes;  thus,  it  can  be  computed  once  for  the  entire  data  set.  As  a  result,  there  is  not  an 
appreciable  computational  cost  in  calculating  atesl(x,y)  . 

Expanding  Equation  (28)  and  pulling  the  division  out  of  the  summation  yields 

whigh  zhigh  whigh  zhigh 

X  X  z^hpsF  ( w  +  x,z  +  y )  XX  ,1psf  (w  +  x,z  +  y) 


Korr(X,y)  = 


Xx,y) 


B(x,y)~ 


Xx,y) 


(58) 


Since  the  PSF  has  been  nonnalized  over  the  window,  the  summation  of  the  PSF  over  the 
window  is  identically  one,  so  Equation  (58)  can  be  simplified  to: 


Lcorr(X,y) 


whigh  zhigh 

z  z  d  ( w,  z)hPSF  ( w + x,  z  +  y) 

W=W ’low  Z=Zlow 


° corr(x,y ) 


B(x,y) 

^corr(x,y) 


(59) 


Since  there  are  N2  pixels  in  the  window,  one  operation  on  a  single  pixel  of  the 
summation  requires  N 2  operations,  so  multiplying  d(w,z)  and  h/JS/,  ( vv  +  x,  z  +  _y)  and  then 

performing  the  summation  requires  two  operations  per  pixel,  or  2N2  - 1  operations  in 
total.  Dividing  the  summation  by  the  standard  deviation  and  then  subtracting  the  ratio  of 
the  background  to  the  standard  deviation  requires  only  three  more  operations.  Thus,  the 
additional  computational  cost  of  correlator  over  the  point  detector  converges  to 
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(60) 


(2N1  -l)  +  3  «  2N2 

Recalling  that  the  cost  of  implementing  the  point  detector  is  3V2 ,  the  cost  of  the 
correlator  is 

3N2  +2N2  =5N2  (61) 

Thus,  the  correlator  requires  67%  more  operations  than  the  point  detector  when  N  is 
chosen  in  accordance  with  current  practice  [15]. 

V.  Conclusions 

1.  The  Correlator  Performs  Better 

The  correlator  outperforms  the  point  detector  in  tenns  of  signal  to  noise  ratio  (SNR) 
in  nearly  all  frames.  More  importantly,  the  correlator  yields  a  higher  probability  of 
detection  of  dim  objects  near  the  threshold  of  detection.  Based  on  the  five  low-passed 
probability  of  detection  curves,  the  correlator  averages  a  75%  chance  of  detection  when 
the  point  detector  has  a  50%  chance  of  detection,  and  the  correlator  averages  a  50% 
chance  of  detection  when  the  point  detection  has  a  20%  chance  of  detection.  Thus,  it  is 
possible  that,  without  changes  made  to  the  way  data  is  collected  and  without  an  increase 
in  the  likelihood  of  false  alann,  dim  objects  will  be  detected  more  frequently  and  even 
dimmer  objects  will  be  detectable. 

2.  Minimal  Cost  to  Implement 

Very  few  changes  are  necessary  in  implementing  the  correlator  versus  the  point 
detector.  The  correlator  only  requires  67%  more  calculations  than  the  point  detector. 
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Additionally,  the  additional  calculations  do  not  extend  the  depth  of  the  pipeline  required, 
and  so  the  additional  computational  cost  can  easily  be  accommodated  by  adding 
hardware  in  parallel. 

3.  Limitations  and  Possible  Future  Work 

The  output  of  the  correlator  changes  less  than  the  output  of  the  point  detector  when 
other  sources  of  illumination  are  in  the  window.  Should  a  particular  task  require  this 
sensitivity,  some  adjustments  may  need  to  be  made.  This  also  suggests  that  the  correlator 
may  be  better  able  to  differentiate  between  closely  spaced  objects.  The  detection  of  dim 
objects  in  the  presence  of  other  objects,  especially  brighter  objects,  could  be  a  very 
fruitful  area  of  study. 

Both  methods  still  rely  on  a  Gaussian  model  of  the  emission  or  reflection  of  photons 
from  the  dim  object,  which  is  inaccurate.  In  truth,  photon  arrival  is  best  described  by  the 
Poisson  distribution.  The  noise  of  the  amplifiers  in  the  circuit  does  create  Gaussian 
noise,  but  any  noise  that  results  in  a  voltage  lower  than  signal-ground  is  registered  by  the 
analog-digital  converters  (ADCs)  as  a  zero.  The  Poisson-based  log-likelihood  ratio  test 
detector  that  Capt  Peterson  developed  could  be  adapted  to  collections  from  the  SST  that 
lacks  a  known  background,  but  establishing  a  metric  for  comparing  Poisson-  and 
Gaussian-based  detectors  remains  a  challenge. 

The  method  of  selecting  a  PSF  used  in  this  research  is  labor  intensive  and  highly 
subjective.  Furthennore,  since  the  imaging  plane  is  under-sampled  for  the  typical  seeing 
conditions,  it  is  likely  that  the  perfonnance  of  the  correlator  suffers  from  the  affects  of 
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aliasing.  The  appearance  of  structure  in  the  periodic  dips  of  the  SNR  of  the  detectors 
supports  the  theory  that  aliasing  is  happening  as  small  shifts  of  the  satellite  on  the 
imaging  plane  take  place. 

4.  Summary 

The  correlator  offers  significant  perfonnance  advantages  for  the  LINEAR  and  SST 
programs  for  little  to  no  cost  depending  on  the  available  computing  resources.  In 
particular,  the  correlator  is  able  to  capitalize  on  the  unprecedented  light-gathering  ability 
of  SST  and  enable  it  to  see  hither-unknown  or  undiscoverable  dim  objects. 
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